NASA CONTRACTuk 
REPORT 


NASA CR- 144294 


SPACE PROCESSING FLOAT ZONE THERMAL ANALYSIS 

By J. T. Poison and D. M. Anderson 

The Boeing Company 
Seattle, Washington 


May i976 


(NASA-CI— 144294) SPACE PROCESSING FLCAT 
ZONE THEENAL ANALYSIS (Boeing Co., Seattle, 
Hash.) 62 p HC $4.50 CSCL 20B 


G3/76 


Prepared for 

NASA-GEDRGE C. MARSHALL SPACE FLIGHT CENTER 
Marshall Space Flight Center, Alabama 35Si2 


N7O-25059 

Unclas 

28245 






1 


1 


t 


technica 


1 REPORT HO. 

NASA CR-144294 


2. government ACCESSIC-' no. 


« title and subtitle 

Space Processing Float Zone Thermal Analysis 


.■t.-iiQPT ST^NOARf! TIT I £ PAGE 
3 • C'-'IFNT' catalog KO 1 


5. REPORT DATE 

May 1976 

A,. PCRTCRMING ORGANIZATION C(>OE 


7. AL'THOR(S) 

J. T. Pogson and D. M. Anderson 

9. PERFORMING ORGANIZATION NAME ANO ADDRESS 

The Boeing Company 
Seattle, Washington 


12 SPONSORING AGENCY NAME AND ADDRESS 

National Aeronautics and Space Administration 
Washington. D. C. 2054 G 


t 

J 

j a. pE»FCKMING organization REPriRr e 


|l0. WORK UNIT NO. 

! 

S 1. '■ONTRACT OR GRANT NQ. 


NAS8-31635 

111. t“oe of repor'. 4 p'^‘^'00 covcae: 

Contractor 


1 1. : TORINO AGENCV _c. L 


_■ __ _____ _ _ _ 1 

SUPPLEMENTARY NOTES 
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I ABSTRACT 

j Two thermal analy.sis (BETA) computer program adaptations were prep.Tred to 

! anftlyze phase change histories in crystal specimens; the first program (BETA-CYL) treats 
right circular cylinder configurations and the second, more general, program (BETA-BOB) 
treats a generalized body-of-revolution configuration. A series of computer runs were 
made for silicon material to determine boundary conditions which produce fl.at solidification 
interfaces while, at the same time, minimizing jx'ak temper.Ttures in the molten zone. 

; Flat solidification interfaces arc a goal believed by some investigators to be required to 
’ produce high quality semiconductor materials. The thi‘rmal effects of convection in a molten 
zone were e.xamiried and found to be negligible in comparison to the conduction heat transfer 
of the melt. 
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1.0 INTRODUCTION 

Growth of single crystals in vertical arc-imaging furnaces has been 
successfully demonstrated by various researchers. Early experiments 
with these furnaces using a carbon-arc energy source, encountered 
difficulties in maintaining uniform light intensity for required 
periods of time. Development of high-power Xenon arc lamps 
substantially reduced the instability problem experienced with 
carbon-arcs, and also allowed longer-term operation for crystal 
growth. Experiments have shown that one of the preferred 
approaches for growing single crystals in arc-imaging furnaces is 
the pedestal technique. This technique has alternatively been 
labeled a modified C 2 ochralski, or modified float-zone process. The 
principal advantage of the float-zone process over other crystal 
growth techniques is production of higher purity crystals; achieved 
by eliminating the need for a crucible and the resulting crucible 
contamination. The homogeneity of crystals grown by the float-zone 
process, as with other processes, is affected by circulation patterns 
in the molten semi-conductor material. Carruthers and Grasso have 
extensively studied the stability and fluid flow in floating zone' 
in both a simulated zero-gravity environment and on Skylab. Effects 
of convection patterns on homogeneity of solute distribution and 
electrical resistivity have been reviewed by Benson. Similar 
considerations should apply to float-zone crystal growth, and 
indicate that for the Czochralski process, all convection in the 
melt should be avoided. 

Growth of single crystals in space is expected to substantially 
improve the solute homogeneity of crystals, by reducing or eliminat- 
ing convection within the melt region. It should also be possible 
to form larger diameter crystals because the length of the suspended 
liquid zone can be increased in the absence of gravity. 

Maintaining a large-diameter un i forn,- temper ature melt zone will be a 
key problem in designing a crystal growth experiment using a solar 
concentrator. Accordingly, this study has been directed at thern'al 
modeling of the liquid zone to predict temperature gradients for 

t 

various parametric changes. 



2.0 ANALYSIS 

Two generalized lumped parameter nodal math models were formulated 
to represent a right circular cylindrical, BETA-CYL, and a body of 
revolution configuration, BETA-BOR, with time varying boundary con- 
ditions. The physical system was sectioned into an array of discrete 
volumes and the mass of each volume element was "lumped" at a point 
within the volume which it represented. The paths for heat transfer 
from one mode to another were represented by conductors joining 
the appropriate nodes. The thermal models each contained a 
maximum of 300 nodes. Figure (1) deoicts the right circular cylindri- 
cal nodal math model . 

The thermal math models described above were input into the Boeing 
Engineering Thermal Analyzer Program (BETA). This program solves 
an electrical analog network system using an iterative numerical 
method. Both transient and steady state solutions are treated. A 
logic diagram of the computer program detailing the computational 
steps is given in Figure (2). Several BETA subroutines required 
development to permit numerical simulation of the applied heat flux 
'boundary conditions and the change of phase process. 

The heat pulse subroutine is illustrated in Figure (3). This 
permits the specification of an externally applied heat flux of 
the form: 

„ - „ „ -A|X-Vt -B|^ 

- <'o ' 

where : 

2 

- Peak surface heat flux, Btu/min-in 
X - Distance from reference end of crystal, inches 
V - Pulse velocity away from reference end of crystal, in/min 
t - Time elapsed since start of pulse movement, min 
B - Starting position of center of pulse measured from reference 
end of crystal , inches 
A - Pulse width parameter 
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Values of A of 100 and 4 were used, the relative widths of the heat 
pulses can be seen in Figure 4. 

A similar cooling subroutine was also formulated to represent 
surface cooling. 

During a phase change each node must receive or release a quantity of 
energy equal to its latent heat of fusion. This energy exchange 
accounting system was simulated by a phase change subroutine, PHAS. 
Figure (5) illustrates the logic system. Each node is surveyed to 
see if it is at or undergoing a phase change. If it is undergoing 
a phase change then an accounting system keeps the nodal temperature 
constant at the melting temperature until the energy exchange matches 
the node latent heat of fusion. 

This study did net directly consider the transparency of silicon 
to thermal radiation; however, the "no convection" and "Infinite 
convection" extremes considered for energy transport in the moltf’n 
zone could be expected to bound the effects which might be observed 
had internal radiation been considered. Furthermore, with a 
cylinc'^ical crystal, any radiation interchange between the 
surroundings and crystal internal regions due to transparency can 
only tend to reduce radial temperature gradients and consequently 
give a more uniform axial flow of thermal energy which is required 
to obtain flat interfaces. Radiation interchange in the axial 
direction, between internal crystal regions, was not expected to 
be very significant as temperature differences within the crystal 
are much less than temperature differences which exist between 
the crystal and its surroundings. 
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3.0 DISeUSSION OF RESULTS 

A total of 82 runs were made to de-bug, checkout, and perform 
production analyses for this study. Only the runs shown in Table 1 
contained significant results. The other runs, although more 
numerous, were less significant in terms of expended machine time. 

The first block on this table describes the crystal configuration 
for each run and refers to a right circular cylinder unless other- 
wise stated. The right circular cylinder crystals were considered 
to consist of 30 wafers of equal thickness as defined by planes 
normal to the cylinder axis (see "Wafer Thickness" column). 

The surface heat flux refers to the peak surface heat flux of a 
Gaussian shaped pulse when the word "Gauss" appears in the "Pulse 
Width" column. If the word "Uniform" appears in the "Pulse Width" 
column, the surface heat flux value applies over the entire heated 
length, and the number appearing in the "Pulse Width" column 
is the heated length, in inches, of a section of cylinder centered 
along the cylinder axis. 

For some runs, the peak surface heat flux, q^, was reduced linearly 
with time; the q^ limits where applicable, are described in the 
"Surface Heat Flux" column. In any case, when the run time,»i.e. 
"Real Time" in the "Time Data" block, exceeds the time shown in the 
"Heat Time" column, the heat is reduced to zero after the "Heat 
Time" shown. 

In the "Number of Nodes" column, "300" indicates that each wafer was 
considered to consist of 10 rings, the temperatures of which were 
computed as node temperatures by the program. The number "120" 
indicates each wafer was considered to consist of 4 rings. 

Under the "Cooling Parameters" block, the "Sink Temperature" and 
"Script F" define radiation rates from the redial surface of the 
crystal for the axial length of surface specified in the "Cooled 
Area" column. The ends; normal to the cylinder axis, were con- 
sidered insulated for all runs. 
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Except for Runs 35, 72, 75, 80, 81, and 82, the initial temperature, 
uniform throughout the crystal, was one degree Fahrenheit below 
the melting point. This was done to prevent wasted computer time 
: * raising the crystalto the melting point. Runs 72, 75, 80, 81, and 

82 were intended to determine steady state conditions and estimates 

I ' of the final temperatures were used as initial temperatures for the 

I 

f analysis to conserve machine time, i.e., initial temperatures above 

the melting point were used in the heated reoion, and temperatures 
somewhat below the melting point were used in the cooled 
regi on . 

Run 35 was a check of the program against an analytical solution for 
a cylinder, with uniform constant surface heat flux, using an 
initial uniform temperature of zero; excellent agreement was 
obtai ned. 

3.1 Cylindrical Melt Zones 

The initial part of this study consisted of analyzing cylindrical 
specimens (BETA-CYL program) on a transient basis using a constant 
strength heat source moving in the axial direction with ends 
cooled by radiation to a 70°F environment as shown below: 


Heat 



The source of the incident radiant energy was not modeled since 
it was considered that it could be best treated after determining 
the boundary conditions required to produce flat interfaces. The 
following parameters were varied in the transient analyses: heated 
length, cooled length, specimen diameter, constant cource strength 
(moving heat sources), and transient source strength (moving and 
fixed heat sources). It was determined that an adiabatic section 
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between the heated and cooled sections would be beneficial in pro- 
ducing flat interfaces. Additionally, it was found that excessive 
machine time was required with a transient-type analysis due to the 
rather low *‘,ol id i f i cati on rate required for growth of high-perfection 
crystals («1 cm/hr). A preferable analytical technique, in com- 
parison to the transient analysis, for treating the continuous zone 
refinement/growth along a rod was to obtain the steady-stage solutions 
for configurations with stationary boundary conditions which give 
flat interfaces. If an extremely slow feed rate of the specimen is 
then introduced a'ter steady state is achieved, a negligible 
change in the shape and location of the interfaces would occur with 
respect to a fixed reference system. Machine time required to 
solve for steady-state conditions is nominal provided that reasonable 
estimates of the final conditions are provided as inputs. 

Figure 6 describes the conditions and results for Run 71. Excellent 
flat-interface solidification conditions were obtained and are shown 
in Figure 7. The physical state of each ring of each wafer has been 
represented in the computer printout by an "0", "1", or "2". Each 
column of figures represents a wafer of the crystal as defined in 
Table 1, The radial width of each wafer, as defined in the program, 
varies so that the volumes of each ring are equal regardless of the 
nominal ring radii. The figures are, therefore, not to scale in 
the radial direction. Since all wafers are the same thickness 
in the axial direction (for the right circular cylinder configuration) 
the printouts are to scale in the axial direction. The axis of the 
cylinder is indicated by a centerline at each end of the crystal. A 
"0" indicates that the state of that node is solid. A "1" indicates 
a molten state, and "2" indicates a partially melted state at the 
melting temperature. A column of zeros adjacent to a column of twos 
or ones indicates a flat solidification interface within a one-wafer- 
thickness tolerance. These printouts were used as a visual tool to 
allow rapid scanning of run results for each case. More exact 
definition of interface shapes can be determined by inspection 
of node freezing fraction data in the computer result numerical 
output . 
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The data printed at the left of each state are: 

Time - minutes from start of analysis 

6. - maximum axial temperature gradient, *F/inch 

''max 


Sb " maximum radial temperature gradient, ®F/inch 

•^MAX 

- maximum temperature, ®F 


A melting point temperature of 2573®F was assumed for silicon. 
Conductivity and specific heat values employed as a function of 
temperatures are shown In supplementary Tables 2 and 3. 


Run 71 was for a 1-inch diameter crystal, 7.5 Inches long with a 

stationary uniformly irradiating heat source applied to the central 

2 Inches of specimen. The source strength decayed linearly from 
2 

20 Btu/mln-ln to 0 In the first 2 minutes of run time. For 
approximately 7 minutes after the heat source was shut off, the 
central molten zone solidified with virtually flat Interfaces. 

The solidification rate was approximately 40 cm/hr, which Is higher 
than current state of the art. The solidification process, 
however, could be Increased to any desired duration by reducing 
the cooling rate at the crystal ends.e.g., through the use of louvers 
or radiation shields. The simulated high solidification rate was 
used here to conserve machine time. A lower cooling rate would be 
accompanied by lower temperature gradients in the crystal compared 
to those shown on Figures 6 and 7. 


Figures 8 and 9 contain the computation for runs 75 and 82 
respectively. A section through one half of the length of crystal 
Is shown with each node defined by a rectangle (to scale), with 
resulting steady state temperatures given at each node in degrees F. 
The melt Interface shape was sketched In based on freezing fraction 
data from the numerical output of the program. 


A nearly converged solution using the steady state analysis technique 
was obtained with the case analyzed In Run 75; the corresponding 
conditions and results arc shown on Figure 8. Although there was 
a slight deviation from flatness at the Interface, It Is believed 
that this could be eliminated by a slight reduction In the rate of 
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source Irradiation which would move the Interface Into the Insulated 
(shielded) region where the thermal flux lines are more nearly 
straight. 

To bound the effects of convection In the molten zone, an "Infinite 
convection" analysis was performed for the same configuration as in 
Run 75. This condition was Implemented by assigning a cup-mixinc 
temperature to the melt at each Iteration which is the effect of 
Instantaneously distributing throughout the melt zone all Incident 
energy at the surface. The results are shown In Run 82, Figure 9. 

A virtually flat Interface was obtained In almost the same position 
as with no convection. By comparing figures 8 and 9 It Is seen 
that there Is little difference In the two solutions except that 
the melt temperature with convection Is uniform at a value somewhat 
lower than the peak value observed with no convection. The degree 
of convergency to a steady state solution can be assessed by noting 
how closely the total cooling rate approaches the constant heating 
rate. For Runs 75 and 82 the errors are 3% and 1" respectively, which 
is sufficient convergence for preliminary purposes. 

Additional transient results. Runs 29, j6, 39, 40, 44, 49, 55, 56, 

60, 62 and 65, similar to those presented for Run 71 are shown In 
Figures 10 through 20. The significant results for each of these 
runs are shown In the "Comments" column of Table 1. 

Crystal specimen dimensions and boundary conditions are shown on 
Figures 4, 21, and 22 for Runs 49, 55, 56, 60, and 62. Of parti- 
cular Interest are the heat pulse shapes and positions on the speci- 
men as a function of time which cannot be readily visualized from 
Table 1 data. 

For Runs 60, 62, and 65 maximum axial and radial temperature 
gradients, maximum temperature, and peak surface heat flux are 
shown as a function of time on Figures 23, 24, and 25 respectively. 
The Irregularities In the temperature and gradient curves are due 
to (1) the results of a lumped parameter thermal analysis, e.g. 
at the Instant when a node melts or solidifies, a relatively 
sudden adjustment In adjacent node temperatures occurs, and (2) 
the location of the peak temperature and gradients moves around 
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in the crystal , whereas, the plots show only the peak values 
regardless of location in the crystal. 

3.? Body of Revolution Melt Zones 

The body-of-revol ution program (BETA-BOP) has the capability to 
analyze configurations where the specimen has a shape more complex 
than a right circular cylinder. Possible configurations of interest 
are shown on Figure 26: (a) vertical cylinder float zone with gravity 
(b) float-zone in Zero-G with rotation, and (c) crystal growth from 
different sized feed stock. An example of case (c) was analyzed 
for 6-inch and 12-inch rods meeting in a molten zone with an assumed 
curved section transition; the configuration and boundary conditions 
are shown on Figure 27. Estimates of steady state conditions used 
for program input are shovm in Figure 23. Conditions after 2.35 
minutes are shown on Figure 29. Steady state convergence is 
approached with a good flat solidification interface on the 12-inch 
end. The interface shapes were sketched in on Figure 29, based 
on melt fraction data from the pregram output for the respective 
interface nodes. 



Figure 27 shows the locations of the nodes for Run 80 as limited 
by the ability of the node size to be resolved into cells on the 
printout page. The dashes and asterisks are used i;, an alternate 
fashion for each wafer to define the node bcundaries for each wafer. 
Each wafer has 10 ring shaped nodes, however on the 6" diameter 
end, not all are shown due to the resolution limitation described 
above. The radial interfaces between the nodes in regions of body 
diameter change are considered by the program to be conical surfaces. 
The uniform irradiation rate shown for the heated area is based 
on actual surface area, not a correct area based on some assumed 
ray direction. The shape shown on the printout is approximately 
to scale. 

Figure 28 shows the initial conditions for Run 80: an asterisk 

indicates a temperature abeve the melting point, a dash indicates 
a temperature below the melting point. The difference between 
this figure as compared to Figure 7 is that the symbols at each node 
correspond to the node state rather than alternating to define node 
boundaries. Final conditions are shown for Run 80 on Figure 29 

0 

original page is 
OP POOE QUALITY 


J L t _ I i 

after the Indicated elapsed time, given the Figure 28 initial 
conditions. Dashes and asterisks are as described for Figure 28 
and an "X" symbol denotes a node partially melted at the melting 
temperature. Melt interface shapes shown are based on freezing 
fraction daU for the melting nodes. Temperatures at selected 
points were shown for reference. Heating and cooling rates were 
shown as a guide to the deviation from steady state, i.e. if a 
true steady state had been reached, heating and cool ing rates 
would be equal. This condition could be achieved by allowing 
the run to proceed in time, at a corresponding cost in machine time. 

Figures 30, 31, and 32 show data similar to that of Figures 27. 28, 
and 29 which apply to Run 72; neither steady state conditions nor 
flat interfaces were obtained. 
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CONCLUSIONS 


Thermal math models of a silicon float zone crystal growth process 
were formulated. Transient temperatures and melt zones were 
determined for a variety of simulated external heat flux conditions. 
It was found that nearly flat interface solidification conditions 
could be simulated by varying the applied heat flux« as in cases 
•65 and 71. Additionally, convective affects were found to be 
negligible in reducing temperature gradients. 
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SILICON THERMAL CONDUCTIVITY 


TABLE 


TEMP - ' V 


Thermal Conductivity Btu/nun-in-°F 


32 

BO. 3 

170.3 
::v0.3 

440.3 
620. 3 

800.3 
980. 3 
1160 
1340 
1520 
1880 
2440 
2573 


.1348 
.1187 
.0955 
.0793 
.0611 
.04 97 
.0408 
.0339 
.0287 
.0250 
.0223 
.0196 
.0178 
.0176 
.0176 
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SILICON SPECIFIC HEAT 


TABLE 3 

TEMP ~ *F Specific Heat - Btu/LB-**F 


0 

.157 

200 

.184 

400 

.197 

600 

.203 

800 

.208 

1000 

.211 

1200 

.216 

1400 

.220 

1600 

.224 

1800 

.228 

2000 

.234 

2200 

.242 

2400 

.245 

2573 

.249 

2600 

.250 

io6 

.250 


250 






BETA PROGRAM - CRYSTAL MELT/FREEZE TRANSIENT STUDY 





FIGURE: 2 
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HEAT PULSE SUBROUTINE (QSFC) 


Q(Z,t,V,A,B), Btu/inin-ln 


iimiiyiyimiiiuiiiiiBi 



PULSE WIDTH a 1/A 


CYLINDER 
SURFACE 
NODES (TYP.) 


FLOW DIAGRAM 


CALL QSFC (Qo.t,V.A,B,Q^) 


50 TO FIRST SURFACE NODE 


COMPUTE HEAT . 
INPUT TO NODE i’ 


*= itDQq y g-A(Z-vt-B)‘ 


GO TO NEXT 
SURFACE NODE 


^ WAS ^ 
THAT THE 
LAST SURFACE 
. NODE ? . 


RETURN TO 
BETA PROGRAM 


COOLING SUBROUTINE ISSIMILIAR 


FIGURE: 3 




LOCAL SURFACE HEAT LOCAL SURFACE HEAT 

FLUX (BTU/MIN-IN 2) FLUX (BTU/MIN-IN .2) 






DISTANCE FROM END OF CRYSTAL (INCHES) 

Figure 4 

RUN 49 & 55 - LOCAL SURFACE HEAT FLUX PROFILES 
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PHASE CHANGE SUBROUTINE (PHAS) 













MAX MEtT TEMPERATURI 
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Figure 8. Run 75 STEADY STATE-NO CONVECTION 
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Figure 9. Run 5? STEADY STATE-INFINITE CONVECTION 
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Figuny 9 . Hun 62 {Continued) 
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LOCAL SURFACE HEAT LOCAL SURFACE HEAT 

FLUX (BTU/r.:iN-IN .2) FLUX (BTU/MIN-IN .2} 


RUN 56 CRYSTAL DIA» 1.5 INCH 
ENDS INSULATED 

TIME - MIN 

i ii 



Figure 2 1 

SURFACE HEAT FLIX RUN 56 & 60 

47 
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(•) VERTICAL CYLINDER WITH GRAVITY 



(b) FLOAT ZONE IN ZERO-G WITH ROTATION 



ic) FLOAT ZONE BETWEEN DIFFERENT DIAMETERS 


FIGURE; 26 . Typical Body-of-Revolution Configurations 
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MAP SHOWING NOO£ LOCATIONS AND BOUNDARY CONDITIONS 



FIGURE: 27 BODY OF REVOLUTION NODE MAP 
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FIGURE: 28. ^unSO - MELT ZONE - TIME - 0 
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CONDITION AT TIME » 2.8S2 MIN 
HEATING RATE - 661 8TU/MIN . 
COOLING RATE - 594 BTU/MIN 
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MAP SHOWING NODE LOCATIONS AND BOUNDARY CONDITIONS 
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INITIAt. CONDITIONS - TIME - 0 
HEATING RATE - 084 BTU/MIN 
COOLING RATE - 1013 BTU/MIN 



Figure MELT ZONE TIME • 0 RUN 72 
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